##=============================================================================
## Appendix Figure 19
##=============================================================================

##-----------------
# clear environment
rm(list=ls())
options(stringsAsFactors = FALSE, scipen = 999)
# source("R/functions.R")

seed <- sample.int(.Machine$integer.max, 1)
set.seed(seed)

ipak <- function(pkg){new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}

packages <- c("tidyverse", "ordinal", "broom", "ggeffects", "lmtest",
              "mice", "flexplot")

ipak(packages)

##---------
# Load data
#setwd("/Users/austinknuppe/Library/CloudStorage/Dropbox/Ukraine2022WartimeSurvey/Paper_peace/final_version_oct_2024/replication-scripts")
load("clean_ukraine_data.RData")

##----------------------------------------
# Sensitivity Analysis with Imputed Values

# https://www.gerkovink.com/miceVignettes/Sensitivity_analysis/Sensitivity_analysis.html

##-------------------------------------------------------------------------
# Let's start by assuming MCAR -- DKs and Ref have the same DGP as observed

dat <- dat %>% 
  mutate(save_lives = as.factor(save_lives))

fit <- clm(save_lives ~ survival_mindset, data = dat)
tidy(coeftest(fit, conf.int = TRUE, conf.level = 0.95))
tidy(coefci(fit)) # observed 1.73 [1.47, 1.98]

g <- ggpredict(fit, terms = "survival_mindset")

g <- g %>% 
  as.tibble() %>% 
  dplyr::select(-group) %>% 
  mutate(response.level = case_when(
    response.level == 1 ~ "Much more important to keep \n territory than to save lives",
    response.level == 2 ~ "Somewhat more important to keep \n  territory than to save lives",
    response.level == 3 ~ "Somewhat more important to \n save lives than to keep territory",
    response.level == 4 ~ "Much more important to \n save lives than to keep territory")) %>% 
  mutate(response.level = factor(response.level, 
                                 levels = c("Much more important to keep \n territory than to save lives",
                                            "Somewhat more important to keep \n  territory than to save lives",
                                            "Somewhat more important to \n save lives than to keep territory",
                                            "Much more important to \n save lives than to keep territory")))


ggplot(g, aes(x=x, y=predicted, colour=response.level)) + 
  geom_ribbon(aes(ymin=conf.low, 
                  ymax=conf.high, 
                  fill=response.level), 
              alpha=.15, 
              col="transparent") + 
  geom_line() + 
  theme_bw() + 
  theme(legend.position="bottom") + 
  labs(x="Survival Mindset", 
       y = "Predicted Probability", 
       colour = "Save Lives (Ordinal)", 
       fill = "Save Lives (Ordinal)")

# Perform a dry run (using maxit = 0) in mice. 
#List the number of missing values per variable.
ini <- mice(dat, maxit = 0)
ini$nmis
md.pattern(dat)
fluxplot(dat)

# observed values of IV
mean(dat$survival_mindset, na.rm = TRUE) # 0.32
sd(dat$survival_mindset, na.rm = TRUE) # 0.25
0.32 + qnorm(0.975) * 0.25
0.32 - qnorm(0.975) * 0.25

summary(dat$survival_mindset)
hist(dat$survival_mindset)

##----------------------------------------------
#  Generating imputations under the δ-adjustment
delta <- c(0, -0.05, -0.10, -0.15, -0.20, 0.05, 0.10, 0.15, 0.20)

imp.all <- vector("list", length(delta))

post <- ini$post

# NB replace for-loop with function for better comp. speed!
for (i in 1:length(delta)){
  d <- delta[i]
  cmd <- paste("imp[[j]][,i] <- imp[[j]][,i] +", d)
  post["survival_mindset"] <- cmd
  imp <- mice(dat, post = post, maxit = 5, seed = i, print = FALSE)
  imp.all[[i]] <- imp
}

#bwplot(imp.all[[1]]) 
imp_1_plot <- densityplot(imp.all[[1]], lwd = 3) # delta = 0
imp_1_plot[["xlab"]] <- expression(paste("Survival Mindset (", delta, " = 0)"))

#bwplot(imp.all[[5]]) 
imp_5_plot <- densityplot(imp.all[[5]], lwd = 3) # delta = -0.20
imp_5_plot[["xlab"]] <- expression(paste("Survival Mindset (", delta, " = -0.20)"))

#bwplot(imp.all[[9]]) 
imp_9_plot <- densityplot(imp.all[[9]], lwd = 3) # delta = 0.20
imp_9_plot[["xlab"]] <- expression(paste("Survival Mindset (", delta, " = 0.20)"))

ggpubr::ggarrange(imp_1_plot, imp_5_plot, imp_9_plot,
                  ncol = 1, 
                  nrow = 3,
                  labels = c("Imputation 1: Missing Completely at Random", 
                             "Imputation 5: Skewed -1 SD", 
                             "Imputation 9: Skewed +1 SD"))


fit1 <- with(imp.all[[1]], clm(save_lives ~ survival_mindset))
pool.fit1 <- pool(fit1)
summary(pool.fit1, conf.int = TRUE)

fit2 <- with(imp.all[[2]], clm(save_lives ~ survival_mindset))
pool.fit2 <- pool(fit2)
summary(pool.fit2, conf.int = TRUE)

fit3 <- with(imp.all[[3]], clm(save_lives ~ survival_mindset))
pool.fit3 <- pool(fit3)
summary(pool.fit3, conf.int = TRUE)

fit4 <- with(imp.all[[4]], clm(save_lives ~ survival_mindset))
pool.fit4 <- pool(fit4)
summary(pool.fit4, conf.int = TRUE)

fit5 <- with(imp.all[[5]], clm(save_lives ~ survival_mindset))
pool.fit5 <- pool(fit5)
summary(pool.fit5, conf.int = TRUE)

fit6 <- with(imp.all[[6]], clm(save_lives ~ survival_mindset))
pool.fit6 <- pool(fit6)
summary(pool.fit6, conf.int = TRUE)

fit7 <- with(imp.all[[7]], clm(save_lives ~ survival_mindset))
pool.fit7 <- pool(fit7)
summary(pool.fit7, conf.int = TRUE)

fit8 <- with(imp.all[[8]], clm(save_lives ~ survival_mindset))
pool.fit8 <- pool(fit8)
summary(pool.fit8, conf.int = TRUE)

fit9 <- with(imp.all[[9]], clm(save_lives ~ survival_mindset))
pool.fit9 <- pool(fit9)
summary(pool.fit9, conf.int = TRUE)

imp_dat <- tibble(
  var = c("Observed Value",
          "Delta = -0.05",
          "Delta = -0.10",
          "Delta = -0.15",
          "Delta = -0.20",
          "Delta = 0.05",
          "Delta = 0.10",
          "Delta = 0.15",
          "Delta = 0.20"),
  coef = c(1.8329056, 1.6846146, 1.6907694, 1.5733186, 1.6840792,
           1.6965847, 1.7560410, 1.6965847, 1.7882018),
  lower = c(1.4971250, 1.2775136, 1.2665536, 0.9160249, 1.3309989,
            1.0527183, 1.2806453, 1.0527183, 1.3380200),
  upper = c(2.1686862, 2.0917156, 2.1149851, 2.2306123, 2.0371595,
            2.3404511, 2.2314368, 2.3404511, 2.2383836),
  color = c("black", "gray70", "gray70", "gray70", "gray70",
            "gray70", "gray70", "gray70", "gray70")) 

imp_dat <- imp_dat %>%
  mutate(var = factor(var, levels = c("Observed Value", "Delta = -0.05",
                                      "Delta = -0.10", "Delta = -0.15",
                                      "Delta = -0.20", "Delta = 0.05",
                                      "Delta = 0.10", "Delta = 0.15",
                                      "Delta = 0.20")))

labs <- c(expression(paste(delta, " = 0.0")),
          expression(paste(delta, " = -0.05")),
          expression(paste(delta, " = -0.10")),
          expression(paste(delta, " = -0.15")),
          expression(paste(delta, " = -0.20")),
          expression(paste(delta, " = 0.05")),
          expression(paste(delta, " = 0.10")),
          expression(paste(delta, " = 0.15")),
          expression(paste(delta, " = 0.20")))

ggplot(imp_dat, aes(x = var, y = coef, color = color)) +
  geom_hline(yintercept = 1.0, colour = "red", lty = 1, size = 2) +
  geom_linerange(aes(x = var, 
                     ymin = lower,
                     ymax = upper),
                 lwd = 2) +
  geom_point(aes(x = var, y = coef), size = 4) + 
  scale_color_manual(values = c("black", "gray70", "gray70", "gray70", "gray70",
                                "gray70", "gray70", "gray70", "gray70")) +
  scale_x_discrete(labels = labs) +
  scale_y_continuous(breaks = c(1, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6)) +
  labs(y = "Survival Mindset (Odds Ratio)",
       x = "") +
  coord_flip() +
  theme_bw(base_size = 18) +
  theme(legend.position = "none")

rm(list = ls())
##=============================================================================
## End of File
##=============================================================================